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Abstract. The spectral properties of Lai-^Sr^TiOs, a doped Mott insulator with strong Coulomb corre- 
lations, are calculated with the ab initio computational scheme LDA+DMFT(QMC). It starts from the 
non- interacting electronic band structure as calculated by the local density approximation (LDA), and 
introduces the missing correlations by the dynamical mean-field theory (DMFT), using numerically exact 
quantum Monte-Carlo (QMC) techniques to solve the resulting self-consistent multi-band single-impurity 
problem. The results of the LDA+DMFT(QMC) approach for the photoemission spectra of Lai-^Sr^TiOa 
are in good agreement with experiment and represent a considerable qualitative and quantitative improve- 
ment on standard LDA calculations. 

PACS. 71.27.+a Strongly correlated electron systems; heavy fermions - 74.25.Jb Electronic structure - 
79.60.-i Photoemission and photoelectron spectra 



1 Introduction 

At present, the electronic properties of solids are inves- 
tigated by two essentially separate communities, one us- 
ing model Hamiltonians in conjunction with many-body 
techniques, the other employing density functional the- 
ory (DFT) DFT and its local density approximation 
(LDA) have the advantage of being ab initio approaches 
which do not require empirical parameters as input. In- 
deed, they are highly successful techniques for the calcula- 
tion of the electronic structure of real materials Q . How- 
ever, in practice DFT/LDA is seriously restricted in its 
ability to describe strongly correlated materials where the 
on-site Coulomb interaction is comparable with the band 
width. Here, the model Hamiltonian approach is more gen- 
eral and powerful since there exist systematic theoreti- 
cal techniques to investigate the many-electron problem 
with increasing accuracy. Nevertheless, the uncertainty 
in the choice of the model parameters and the technical 
complexity of the correlation problem itself prevent the 
model Hamiltonian approach from being a flexible or re- 
liable enough tool for studying real materials. The two 
approaches are therefore complementary. In view of the 
individual power of DFT/LDA and the model Hamilto- 
nian approach, respectively, it would be highly desirable 
to be able to combine these techniques, thereby creating 
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an enormous potential for future ab initio investigations of 
all real materials, including, e.g., /-electron systems and 
Mott insulators. One of the first successful attempts in 
this direction was the LDA+U method ||, which com- 
bines LDA with a basically static, i.e., Hartree-Fock-like, 
mean-field approximation for a multi-band Anderson lat- 
tice model (with interacting and non-interacting orbitals) . 
This method proved to be a very useful tool in the study 
of long-range ordered, insulating states of transition met- 
als and rare-earth compounds. However, the paramagnetic 
metallic phase of correlated electron systems such as high- 
temperature superconductors and heavy-fermion systems 
clearly requires a treatment that goes beyond a static 
mean-field approximation and includes dynamical effects, 
e.g., the frequency dependence of the self-energy. 

During the last decade a new many-body approach 
was developed which is especially well-suited for the in- 
vestigation of strongly correlated metals - the dynami- 
cal mean- field theory (DMFT) @,|,§. It becomes exact 
in the limit of high lattice coordination numbers and 
preserves the dynamics of local interactions; hence it rep- 
resents a dynamic mean-field approximation. In this non- 
perturbative approach, the lattice problem is mapped onto 
an effective Anderson impurity model with a hybridiza- 
tion function which has to be determined self-consistently. 
To solve the effective impurity problem one can either 
use approximative techniques such as iterated perturba- 
tion theory (IPT) §§ and the non-crossing approxima- 
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tion (NCA) | ,|l^,|ll), or employ numerical techniques like 
quantum Montc-Carlo simulations (QMC) G^, exact di- 
agonalization (ED) |p^,|6|, and numerical renormalization 
group (NRG) Q. In principle, QMC, ED and NRG are 
exact methods, but require an extrapolation: discretiza- 
tion of the imaginary time At — > (QMC), the number 
of lattice sites of the respective impurity model n s — > oo 
(ED), and the parameter for logarithmic discretization of 
the conducting band A — > 1 (NRG), respectively. 

In principle, the main idea of the LDA+U method || 
(i.e., complementing the LDA band structure by a screened 
Coulomb interaction between localized tight-binding or- 
bitals) can be practicall y ap plied with more refined ap- 
proximation schemes p^,|l6|,|F7,18,19|,^0|. Indeed, a calcu- 
lation scheme supplementing LDA with DMFT to include 
dynamic effects was first formulated by Anisimov et al. ]l5[ 
and was used to calculate the photoemission spectra of 
Lai-^Sr^TiO^, a doped Mott-insulator and strongly cor- 
related paramagnetic metal, in connection with IPT flEf 
and NCA (20). The LDA++ approach by Lichtenstein 
and Katsnclson formulates a very similar strategy 
and was recently applied to investigate correlation effects 
in iron |2l],[22|. Both IPT and NCA are approximative 
methods to solve the effective single-impurity problem in 
the LDA+DMFT scheme. They have the advantage of 
being numerically inexpensive, but their reliability, es- 
pecially in the case of multi-band systems with particle 
densities off half-filling, is in principle uncertain. In this 
situation, it is clearly desirable to employ a controlled 
computation scheme to obtain numerically exact results 
from LDA+DMFT. The QMC method is such a scheme 
and was already applied by Lichtenstein and collabora- 
tors to calculate the magnetic excitation spectrum of fer- 
romagnetic iron and the photoemission spectra of 
Sr2RuC>4 23 1 . In the present paper we report on our im- 



plementation of the LDA+DMFT(QMC) technique and 
the results obtained with it for the photoemission spectra 
of Lai-^Sr^TiOs. 



2 Computational scheme 

2.1 Local density approximation (LDA) 

The main problem in combining LDA and model Hamil- 
tonian approaches comes from the fact that their founda- 
tions are very different. The LDA employs a functional of 
the electron density, while Hubbard and Anderson models 
are formulated in terms of localized, site-centered, atomic- 
like orbitals. In order to merge the approaches it is nec- 
essary to write LDA equations in the basis of such or- 
bitals. There is one variant of the LDA method, based on 
tight-binding linearized muffin-tin orbitals in the orthog- 
onal approximation (TBLMTO) which is naturally 
realized in such a basis. The corresponding Hamiltonian 
can be written as 

•f^LDA = (fiilm,jl'm' £ ilm?l'il m 
ilm,jl' m' ,a 

~\~tilm,jl'm'C i il n Cj l r ml ), (1) 



where i, j are site indices, I, m, I', m! are orbital indices, a 
is a spin index, and operators carry a hat. 

Taking this expression as the non-interacting part of 
a multi-band periodic Anderson model, one may comple- 
ment it by a correlation term describing the local contri- 
butions to the Coulomb interaction 
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The prime on the sum indicates that at least two of the 
indices on different operators have to be different. Here, 
Umm' denotes the direct Coulomb integral. The much smaller 
exchange integral and other local contributions of the Coulomb 
interaction have been neglected. Furthermore, non-local 
Coulomb contributions are not considered in H colr . Note, 
that the largest non-local contributions is the nearest- 
neighbor density-density interaction which, to leading or- 
der in Z (Z: number of nearest-neighbor sites), yields only 
the Hartree term [25 which is already taken into account 
in the LDA. In an actual calculation, it is not possible to 
include the local two-particle interaction terms between all 
orbitals appearing in (||) since the number of states grows 
exponentially with increasing orbital quantum number. 
Thus, one usually concentrates on a certain subset of cor- 
relations and treats the influence of the remaining states 
via the use of effective, screened interaction parameters for 
the shells under consideration. In this spirit, we will as- 
sume in the following that it is only necessary to take the 
Coulomb interaction for the d-shell of the transition metal 
ions (i = id and I = Id) explicitly into account; therefore 
the indices il will be omitted. The correlation part of the 
Hamiltonian then acts only on the d-wave functions of 
the transition metal ions. All other valence orbitals will 
be treated as bands of itinerant electrons which are well 
described by the LDA. 

One must take into account, however, that the Coulomb 
interaction is already present in LDA in some averaged 
way. Hence, to avoid double-counting one needs to sub- 
tract this term from the LDA-Hamiltonian. Unfortunately, 
there exists no direct microscopic or diagrammatic link be- 
tween the Hubbard model approach and LDA, and it is 
thus not possible to express the LDA energy rigorously 
via the d — d Coulomb interaction parameter U (except 
for the atomic limit where one can make a connection be- 
tween the Coulomb parameter U of the Hubbard model 
and the second derivative of the atomic total energy as a 
function of the number of electrons). While it is known 
that LDA eigenvalues are rather bad approximations for 
excitation energies of systems with strong Coulomb inter- 
actions, the LDA total energy as a function of the number 
of electrons is a much better approximation for the exact 
functional. Furthermore, the values of U obtained from 
LDA calculations often agree well with experimental data 
and more rigorous calculations. 

Therefore one may expect that the LDA part of the 
Coulomb interaction energy is well approximated by the 
averaged value of the Coulomb interaction energy in (2) 
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Here, U is the mean value of the Coulomb interaction 
and may be obtained from a first-principles constrained 
LDA calculation [p6| or from experiment, e.g. high-energy 
spectroscopy; n d is the total number of d-electrons. 

In LDA, one-electron energies are defined as deriva- 
tives of the total energy as a function of the occupation 
numbers for the corresponding states. Hence, the one- 
electron energy level for the non-interacting d-states is 
obtained by [pL 5[] 

s° d := ^(^lda - E COII ) = ef A U(n d - ±) (4) 

with e^ DA := ^^--Elda, and £"lda being the total energy 
calculated from TJlda- 

Then, the new non-interacting Hamiltonian will have 
the form 

-^LDA = ^ ] (3ilrnjl'm'£il m ™ilm ^ tilrnJl'm'Cii m Cjl' m ')i 
ilm ,jl' m' ,a 

where e\ lm = e am (1 - £z,j d ) + e d 5i^ d . In reciprocal space, 
the matrix elements of the operator Hl D a are given by: 

(#LDA O^gim^'Z'm' = (#LDA (k))gim,g'Z'm' 

-S q l m .,q'l'm'S q l.q d l d U(n d ~ -). (5) 

Here q is an index of the atom in the elementary unit 
cell, (-ffLDA(k)) g ; mi9 'i' m ' is the matrix element in k-space 
of -£/lda, and q d denotes the d-atoms in the unit cell. 
This non-interacting part ff{| DA , taken together with the 
interaction part (|J), 

H = i?LDA + -^corr (6) 

forms the ab initio Hamiltonian H for a particular mate- 
rial under investigation. 

2.2 Dynamical mean-field theory (DMFT) 

In general, the investigation of the correlated-electron Hamil- 
tonian H is too complicated to allow for an exact solution 
or even a numerical investigation with more than about 
10 sites. Here, the DMFT is a powerful approxima- 

tion scheme which takes into account electronic correla- 
tions and, in particular, correctly describes the formation 
of a coherent quasiparticle band and incoherent Hubbard 
bands. 

The DMFT maps the lattice problem onto a single- 
site problem, which is equivalent to a multi-band single- 
impurity Anderson model, with the self-consistency con- 
dition [jl5| (the fc-integrated Dyson equation) 

Gqlm,q' I'm' i z ) = / ~7T~ \zO~ q lm, q 'l'm' (-^LDA \~) )qlm.q' I'm' 
J Vb 

-Sql, qi u£qlm,q'l'm'{z)]~~ ■ (7) 

Here, [•••] implies the inversion of the matrix with el- 
ements n (—qlm), n'(=q'l'm'), and integration extends 



over the Brillouin zone with volume Vb- In the present 
study we consider a cubic-crystal structure and assume 
the tig orbitals to be interacting. Due to the high sym- 
metry of the crystal, these three t^g orbitals are degener- 
ate. Without symmetry breaking, the Green function and 
the self-energy remain degenerate, i.e., G q i mA n' m '(z) = 

G(z)8qlm,q'l'm' and Sqlm,q'l'm'{z) = £{z)8 q l m<q H' m i for I — 

l d and q = q d (where l d and q d denote the Ti ti g orbitals). 
If the partially filled band under consideration is well sep- 
arated from other bands, as in the case of Lai-^Sr^TiOa 
and other transition metal oxides, one can describe the 
physics of the partially filled band at the Fermi energy ap- 
proximately by an effective three band Hamiltonian ff^Mi 
e.g., by downfolding to a basis with ti g orbitals only. One 
obtains (indices I = l d and q = q d suppressed): 

— [(z - E(z))6 m , m > - (^daWW]- 1 . (8) 

Due to the diagonal structure of the self-energy the in- 
teracting Green function can be expressed via the non- 
interacting Green function G°(z): 

G(z) = G°(z-£(z)) = [du> Q£ , (9) 

J Z — Zj (z) — U) 

Thus, it is possible to use the Hilbert transformation of the 
unperturbed LDA-calculated density of states p°(ui) and 
we do so in the following. This approximation is justified 
if the hybridization between the ti g orbitals and the other 
orbitals is rather weak as in the case of LaTi03. 

The DMFT single-site problem depends on = 
G _1 + U and can be formulated in terms of Grassmann 
variables if> and ip*. For the local Green function at a 
Matsubara frequency u v = (2v + 1)tt//3, orbital index to, 
and spin a one obtains: 

GZm = -~J DMPlfltVCe- 4 !*^*' 5 " 1 ! (10) 

where the single-site action A is given by 

A[i/>,r,G- 1 ]= £ ^ImiQZmY^lm 
i>,cr,m 

~\ E / dTrm(rWm(rWm'(rWm'(r) .(11) 
(mo-)^(mcr'Jo 

2.3 Quantum Monte-Carlo method (QMC) 

To solve the effective single-site problem defined above we 
will apply the QMC method which allows for a numeri- 
cally exact solution |L2| . To this end, the imaginary time 
is first discretized into A steps of size At (At — 0.25 eV _1 
throughout this work), and r t := tAr with integer t. Then, 
in a second step, the Hubbard-Stratonovich transforma- 
tion 

- ^ eyi v{ Xs ma;m'*'(lPmtlPmt ~ VC' *t VC''t ) }( 12 ) 
s* , ,=fcl 



4 



LA. Nekrasov et al.: Calculation of photoemission spectra of Lai-ajSr^TiOa using LDA+DMFT(QMC) 



is employed, which replaces the interacting system by a 
sum of 2^ 2M " M ) A different non-interacting systems, where 
M is the number of interacting orbitals. Each non- interacting 
system can be solved by a Gauss integration which yields 
a contribution (M s ) _1 det(M s ) to the Green function in 
imaginary time, parameterized by (t, t'), where 



It* = At 2 



-6u> £ A 



ma-.m'a' @ma\rn' a' 



(13) 



Here, a m a:m'a' — 20 (a' — a + S aa > [m! — m]) — 1 has a dif- 
ferent sign if (ma) and (m'a') are exchanged. Since a full 
summation over all non-interacting systems is computa- 
tionally impossible, the Monte-Carlo method is employed 
for importance sampling. Details of the one-band QMC al- 
gorithm in the context of DMFT can be found in Refs. jl2], 
H and for the band-degenerate case in Refs. pTj. 



3 Results and Discussion 

The stoichiometric compound LaTi03 has a perovskite 
crystal structure with a small orthorhombic distortion 
(LTi-O-Ti w 155°) ||. Below T N = 125 K [|§ it 
is an antiferromagnetic insulator |}(| with a Ti magnetic 
moment of 0.45 iib and a small energy gap of approxi- 
mately 0.2 eV flJlf . At doping x = 0.05, Lai_ x Sr x Ti0 3 
undergoes an insulator-to-metal transition and becomes a 
correlated paramagnetic metal with a strongly enhanced 
susceptibility and electronic specific heat coefficient [Q . 

The standard LDA calculation for undoped LaTi03 
yields a density of states (Fig. |l|) which is typical for early 
transition metal oxides with a completely filled 2p oxygen 
band ranging from —8.2 eV to —4.0 eV and a partially 
filled Ti-3o? band above it. Since the Ti-ion has an octahe- 
dral coordination of oxygen ions, the Ti-3d band is split 
into three degenerate ti g and two degenerate e g subbands 
which do not mix. Titanium is three- valent in LaTiOa and 
the corresponding formal ionic configuration is dr. This 
implies a partially filled t^ g subband containing one elec- 
tron with a total capacity of six electrons. The e g subband 
is empty and situated just above the ti g subband. In order 
to simplify our calculations the real orthorhombic crystal 
structure was replaced by a cubic structure with the same 
volume. This approximation leads to a slight overestima- 
tion of the effective bandwidth. 

The tig band is the only partially filled band and is 
well separated from other bands. Therefore in our DMFT 
calculation we took only Coulomb interactions between 
electrons in ti g orbitals into account. The resulting prob- 
lem is equivalent to a Hubbard model with three degen- 
erate bands. The DMFT chemical potential was adjusted 
to yield the doping x — 0.06, i.e. 0.94 electrons in the t 2g 
band. We employ Eq. (||) within constrained LDA j26j, 
i.e., changing the number of f2g-electrons only, to cal- 
culate U. Our LMTO-ASA calculation (TB-LMTO-ASA 
code of Andersen and coworkers [Q version 47) yields a 
fully-screened Coulomb interaction of U — 4.2 eV within 
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Fig. 1. Densities of states of LaTiOa calculated with LDA- 
LMTO. Upper figure: total density of states; lower figure: par- 
tial t2 g (solid lines) and e 9 (dashed lines) densities of states. 



the basis Ti(4s,4p,3d) La(6s,6p,5d) 0(2s,2p) at a Wigner 
Seitz radius of 2.37 a.u. for Ti. Our result has to be com- 
pared to that of Solovyev et al. |54j who obtained U — 3.2 
eV employing the ASA-LMTO method within orthogonal 
representation. This shows that the ab-initio calculation 
of U, which is the interaction between particular fo g ) or- 
bitals, is rather sensitive to the orthogonality of the wave 
functions and, also, to the choice of the orbitals. Unless 
specified otherwise, we will, thus, use U = 4 eV and keep 
in mind that the inherent uncertainty is about 0.5 eV. 

In Fig. |^, the spectral function obtained from our 
LDA+DMFT(QMC) calculation at temperature T ps 1000 K 
is compared with the non-interacting ti g density of states. 
One can see the typical features of the spectra of strongly 
correlated systems: a lower Hubbard band, a well pro- 
nounced quasi-particle peak, and an upper Hubbard band. 
While for the non-interacting case 100% of the spectral 
weight is located in the quasi-particle band, the LDA+DMFT 
spectra are characterized by a spectral weight transfer 
from the quasiparticle band to the Hubbard bands and 
a narrowing of the quasiparticle band. 

The QMC simulations performed in this paper to solve 
the effective multi-band single-impurity model provide a 
numerically exact solution, but require a large computa- 
tional effort which restricted our calculations to temper- 
atures of the order of T = 0.1 eV (« 1000 K). Since this 
technique yields results only at imaginary 
(Matsubara-) frequencies the calculation of the spectral 
function requires an analytic continuation of the spec- 
tral function via, e.g., the maximum entropy method [ p5| . 
Previous LDA+DMFT investigations of the photoemis- 
sion spectra of L&-\ ^Sr^TiQq used a variant of IPT for 
doped multi-band systems ]l5| and NCA Q to solve the 
DMFT-equations, i.e., approximate techniques. In Fig. |^, 
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Fig. 2. Partial t^g densities of states of LaTi03 calculated with 
LDA+DMFT(QMC) (solid lines) and LDA (dashed lines). 

we compare the results obtained within these approxi- 
mations with the numerically exact QMC simulation, all 
at T » 1000 K. One notes that within IPT the shape of 
the upper Hubbard band is not correct. Moreover, there 
is no quasiparticle peak at 1000 K, the reason being that 
IPT underestimates the Kondo temperature considerably 
such that the very narrow quasiparticle peak found at low 
temperatures (see right inset of Fig. ||) disappears already 
at about 250 K. A similarly narrow IPT quasiparticle 
peak was found in a three-band model study with Bethe- 
DOS by Kajueter and Kotliar ||. While NCA comes off 
much better than IPT it still underestimates the width 
of the quasiparticle peak by a factor of two. Furthermore, 
the position of the quasiparticle peak is too close to the 
lower Hubbard band. In the left inset of Fig. ||, the be- 
havior at the Fermi level is shown. At the Fermi level, 
the NCA yields a spectral function which is almost by 
a factor two too small. The shortcomings of the NCA- 
results appear to result from the well-known problems 
which this approximation scheme encounters already in 
the single-impurity Anderson model at low temperatures 
and/or low frequencies |37],|38]]. Similarly, the deficiencies 
of the IPT-results are not entirely surprising in view of 
the semi-phenomenological nature of this approximation, 
especially for a system off half filling. This comparison 
shows that the choice of the method used to solve the 
DMFT equation is indeed important. 

Photoemission spectroscopy of the early transition metal 
oxides provides a direct tool for the study of the elec- 
tronic structure of strongly correlated materials. A com- 
parison of the experimental photoemission spectra JXo| | 
with results obtained from LDA and LDA+DMFT(QMC) 
at 1000 K (4l) are shown in Fig. ||. To take into account the 
uncertainty in U, we present results for U = 3.2, 4.25 and 
5 eV. All spectra are multiplied with the Fermi step func- 
tion and Gaussian-broadened with a broadening parame- 
ter of 0.3 eV to simulate the experimental resolution |{40|| . 
The LDA band structure calculation clearly fails to re- 
produce the broad band observed in the experiment at 
1-2 eV below the Fermi energy Taking the correla- 
tions between the electrons into account, this lower band 
is easily identified as the lower Hubbard band whose spec- 
tral weight originates from the quasiparticle band at the 
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Fig. 3. Comparison of the spectral densities of Lai-^SrajTiOs 
(x — 0.06) as calculated by LDA+DMFT using the approxima- 
tions IPT and NCA, with the numerically exact QMC-result 
at T = 0.1 eV, i.e., approximately 1000 K, and U — 4 eV. Inset 
left: Behavior at the Fermi level including the LDA DOS. Inset 
right: NCA and IPT spectra for the temperature 80 K. 



Fermi energy and increases with U. The best agreement 
with experiment concerning the relative intensities of the 
Hubbard band and the quasi-particle peak and, also, the 
position of the Hubbard band is found for U = 5 eV [[l2| 
p3{ . The value U = 5 eV is still compatible with the ab- 
initio calculation of this parameter. One should also note 
that the photoemission experiment is sensitive to surface 
properties. Due to the reduced coordination number at 
the surface, the bandwidth is likely to be smaller and the 
Coulomb interaction to be less screened, i.e., larger. Both 
effects make the system more correlated and, thus, might 
also explain why better agreement is found for U = 5 eV. 
Besides, the polycrvstalline nature of the sample and, also, 
spin and orbital [f!4| fluctuation, not taken into account 
in the LDA+DMFT approach, could further reduce the 
quasiparticle weight. 

In conclusion, the LDA+DMFT(QMC) approach is 
shown to be a workable computational scheme which merges 
the conventional band structure approach with a recently 
developed many-body technique in combination with a 
numerically reliable evaluation method. Thereby, it pro- 
vides a powerful tool for future ab initio investigations 
of real materials with strong electronic correlations. The 
LDA+DMFT(QMC) approach not only explains the ex- 
istence of the lower Hubbard band in doped LaTi03, but 
also, in contrast to LDA, reproduces the qualitative pic- 
ture of the spectral weight transfer from the quasiparti- 
cle band to the lower Hubbard band, the position of the 
lower Hubbard band, and the narrowing of the quasipar- 
ticle band. 



6 



LA. Nekrasov et al.: Calculation of photoemission spectra of Lai-^Sr^TiOa using LDA-fDMFT(QMC) 



QMC, U=3.2 eV 
QMC, U=4.25 eV 
QMC, U=5.0 eV 
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Fig. 4. Comparison of the experimental photoemission spec- 
trum H, the LDA result, and the LDA+DMFT(QMC) calcu- 
lation for LaTi03 with 6% hole doping and different Coulomb 
interaction U = 3.2, 4.25, and 5eV. 
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